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The existence of a discontinuity in the inverse dielectric constant of the two-dimensional Coulomb 
gas is demonstrated on purely numerical grounds. This is done by expanding the free energy in 
an applied twist and performing a finite-size scaling analysis of the coefficients of higher-order 
terms. The phase transition, driven by unbinding of dipoles, corresponds to the Kosterlitz-Thouless 
transition in the 2D XY model. The method developed is also used for investigating the possibility 
of a Kosterlitz-Thouless phase transition in a three-dimensional system of point charges interacting 
with a logarithmic pair-potential, a system related to efi'ective theories of low-dimensional strongly 
correlated systems. We also contrast the finite-size scaling of the fluctuations of the dipole moments 
of the two-dimensional Coulomb gas and the three-dimensional logarithmic system to those of the 
three-dimensional Coulomb gas. 
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I. 



INTRODUCTION 



Compact U{1) gauge fields in three dimensions are of 
great interest in condensed matter theory, as they arise in 
effective theories of strongly correlated two-dimensional 
systems at zero temperature iSiii Lightly doped Mott- 
Hubbard insulators, such as high-Tc cuprates, are ex- 
amples of systems possibly described by such theories, 
where the compact gauge field emerges from strong lo- 
cal constraints on the electron dynamics iSiSiSii High- Tc 
cuprates appear to fall outside the Landau Fermi liq- 
uid paradigm, and a so-called confinement-deconfincmcnt 
transition in the gauge theories may be associated with 
breakdown of Fermi-liquid and quasiparticles in 2D at 
T = OiSiL^ Obliteration of electron-like quasiparticles 
and spin-charge separation in the presence of interac- 
tions is well known to occur in one spatial dimension. 
However, the mechanism operative in that case, namely 
singular forward scattering, is unlikely to be operative in 
higher dimensions due to the much less restrictive kine- 
matics at the Fermi surfaceiS Proliferation of instantons 
of emergent gauge fields show more promise as a viable 
candidate mechanism. This line of pursuit has recently 
been reinvigorated in the context of understanding the 
physics of lightly doped Mott-Hubbard insulators and 
unconventional insulating states, 

The compact nature of a constraining gauge field on 
a lattice model introduces topological defects defined by 
surfaces where the field jumps by 2tt, forming a gas of 
instantons (or "monopoles") in 2 -I- 1 dimensions. Con- 
sidering the gauge sector only, the interactions between 
these instantonic defects are the same as between charges 
in a 3D Coulomb gas, i.e. 1/r-interactions. Such a 
gas is always in a metallic or plasma phase with a fi- 
nite screening lengthjiiii^ and there is no phase transi- 
tion between a metallic regime and an insulating regime. 
However, in models where compact gauge fields are cou- 
pled to matter fields, the interaction between the mag- 
netic monopoles may be modified by the emergence of 



an anomalous scaling dimension of the gauge field due 
to critical matter-field fiuctuations.^'^ This is the case for 
the compact abelian Higgs model with matter fields in 
the fundamental representation^ 

In Refs. [T3. it was shown that the introduction 
of a matter field with the fundamental charge leads 
to an anomalous scaling dimension in the gauge field 
propagator—. The effect is to alter the interaction poten- 
tial between the magnetic monopoles from 1/r to — Inr. 
The existence of a confinement-deconfinement transition 
in the gauge theory is thus related to whether a phase 
transition occurs in a 3D gas of point charges with log- 
arithmic interactions. However, one should note that 
the legitimacy of a monopole action based on just pair- 
wise interactions has been questioned, particularly when 
viewed as an effective description of an effective gauge 
theory of strongly interacting systems The 3D loga- 
rithmic plasma is however of considerable interest in its 
own right. 

In two dimensions, where — Inr is the Coulomb po- 
tential, it is known that the logarithmic gas experiences 
a phase transition from a low-temperature insulating 
phase consisting of dipoles to a high-temperature metallic 
phase. This is nothing but the Coulomb-gas representa- 
tion of the Kosterlitz-Thouless transition in the 2D XY 
model. In a 3D logarithmic gas, the existence of a phase 
transition is still subject to debateii^iiSiii Renormaliza- 
tion group arguments have been usedi^ to demonstrate 
that a transition may occur, driven by the unbinding of 
dipoles. Others have claimed that the 3D logarithmic gas 
is always in the metallic phase^i^ In a recent paper^ 
large scale Monte Carlo simulations indicated that two 
distinct phases of the 3D-log gas exists; a low-T regime 
where the dipole moment does not scale with system size 
and a high-T regime where the dipole moment is system 
size dependent. Those results do however not determine 
the character of the phase transitions. That will be the 
main subject of this paper. 

The Kosterlitz-Thouless transition in the 2D XY model 
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is characterized by the universal jump to zero of the hehc- 
ity modulusiiS In the corresponding 2D Coulomb gas, it 
is the inverse of the macroscopic dielectric constant e that 
experiences a jump to zero when going from the insulat- 
ing to the metallic phase. According to Ref. such a 
universal discontinuity should also take place for e^^ in 
the 3D logarithmic gas associated with the confinement- 
deconfinement transition. Proving that such discontinu- 
ities exist numerically is a subtle task. The discontinuous 
character of the helicity modulus in the 2D XY model is 
very hard to see in a convincing manner by computing 
the helicity modulus, due to severe finite-size effects. It 
was only recently proven on purely numerical grounds 
that such a discontinuity exists^li in a simple, but yet 
clever manner. By imposing a twist across the system 
and expanding the free energy in this twist to the fourth 
order, a stability argument was used to show that the 
second order term in the expansion, the helicity modu- 
lus, must be nonzero at T^. The proof relies on the ability 
to conclude that the fourth order term is negative in the 
thermodynamic limit, from which the discontinuity fol- 
lows immediately. In this paper, we will repeat this pro- 
cedure, but now in the language of the 2D Coulomb gas. 
In addition to confirming the results of Minnhagen and 
Kim, the method which we develop here could be suit- 
able for proving the possibly discontinuous behaviour of 
in the 3D logarithmic gas. This is a main motivation 
for translating the procedure of Ref. |23 to the vortex 
language, since the 3D logarithmic gas is not the dual 
theory of any simple spin model. After having demon- 
strated the discontinuity in the 2D Coulomb gas, we go 
on to apply the method on the 3D logarithmic gas. We 
also compare the scaling with system size of the mean 
square dipole moment for these logarithmic plasmas, and 
contrast the results with those of the 3D Coulomb gas. 
This is important, since the mean square dipole moment 
does not scale with system size below a certain temper- 
ature for the logarithmic plasmas.^® This indicates that 
two phases exist, where the low-temperature regime con- 
sists of tightly bound pairs. However, the results for the 
3D Coulomb gas are qualitatively different, in accordance 
with the fact that such a low-temperature phase is absent 
in that case. 



II. 



MODEL 



The Hamiltonian of the 2D XY model on a square 
lattice modified with a twist T{x,y) is 



H 



XY 



— — J ^ cos( 



27rry -T), 



(1) 



where is the displacement between the nearest neigh- 
bour pairs to be summed over. We set the coupling con- 
stant J to unity. The volume of the system, i.e. the num- 
ber of lattice points, is L^, and the angle 6i is subject to 
periodic boundary conditions. In the Villain approxima- 
tion, a duality transformation leads to the Hamiltonian 



(2) 



where rrii are point charges on the dual lattice, corre- 
sponding to vortex excitations in the XY model. 
is a lattice derivative and e^'^ is the completely anti- 
symmetric symbol. The potential Vij is given by 



q 



iq-(ri- 



cosqx 



(3) 



which has a logarithmic long-range behaviour. Details 
of the dualization are found in appendix ^ As is 
well known, eq. l(2Jl at zero twist describes the two- 
dimensional Coulomb gas (2D CG). In this representa- 
tion, the Kosterlitz-Thouless phase transition of the 2D 
XY model is recognized by a discontinuous jump to zero 
of the inverse macroscopic dielectric constant e^^ at T^. 
We note that the curl of the twist T acts as a modifica- 
tion of the charge field in the 2D CG. 

The free energy of the system is F = —T\nZ, where 
the partition function is given by summing the Boltz- 
mann factor over all charge configurations: 



-H/T 



(4) 



{m} 



Let us write the Hamiltonian in Fourier representation, 
q 

(5) 

where the discrete Fourier transform is defined as in ap- 
pendix IH and A'^e^'i'" = e±''5''Q^q. 



III. 



STABILITY ARGUMENT 



From (HJ, it is clear that F(T) > F{0) in the low- 
temperature phase, i.e. the free energy is minimal for 
zero twist. This inequality is also valid at the critical 
temperature Tc, since the free energy must be a contin- 
uous function of temperature. As a consequence, the 
Taylor expansion 



F{T)^F{0) = J2J2 



OF 



qi 



dT!? 



EE 



qi 
d^F 



T=0 



rpO! 

qi 



a,l3 qiq2 '^-'qi'^^q2 



qi q2 



(6) 



T=0 



can not be negative for any T < Tc. Expressions for the 
derivatives of the free energy with respect to a general 
twist are found in appendix^ Only terms of even order 
will contribute to the series, since mi may take equally 
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many positive and negative values. We are free to choose 
the twist to be 



IV. 



SIMULATION RESULTS 



. A . /27ry 
T(x, y) — — sm — — 



(7) 



where A is an arbitrarily small constant and ?/ = 1 for 
the two-dimensional Coulomb gas. To the fourth order, 
this long-wavelength twist turns © into 



F{T) - F{Q) - — Ck 1 



{CM' 
32 £4^3 



^(mkm_k)^ - ^ (("ikTO-k)^)^ , 



(8) 



where k — {0,2tt/L) and Ck = QkQ^k^k- We recognize 
the paranthesis in the second order term as the dielectric 
response function e^^(k), where k is now the smallest 
nonzero wave vector in a finite system. Note that the 
prefactors in both terms are independent of system size 
as L ^ oo. The crucial argument to use is the same as 
in Ref. If the fourth order term approaches a finite 
negative value at Tc in the limit L — s- oo, the second or- 
der term, e~^(k 0), must be positive to satisfy the 
inequality F{T) > F{0). Furthermore, since we know 
that the inverse dielectric constant is zero in the high- 
temperature phase, it necessarily experiences a disconti- 
nuity at Tc- As we shall see, Monte Carlo simulations 
show that the fourth order term is indeed negative at Tc 
in the thermodynamic limit. 

The argument described above will also apply to a 
three-dimensional gas of point charges interacting via a 
pair potential of some sort, as long as the twist raises 
the free energy in the low-temperature regime. Since the 
curl of the twist T is a vector in that case, one may for 
instance choose the z-component of this vector as the 
perturbing charge in eq. The two three-dimensional 
systems we will consider are the logarithmic gas and the 
Coulomb gas. The expansion © is valid for any sys- 
tem size L. However, to make the change in free energy 
nondivergent as L —^ oo, the twist must be chosen such 
that the terms in the expansion are independent of sys- 
tem size. This is obtained by choosing rj = 2 for the 
logarithmic gas and rj = 3/2 for the Coulomb gas. i] is 
defined in {Tj). In both cases, the second order term will 
be proportional to 



The fourth order term will be proportional to 



e4(k) = ^ 



[rrikm^kf - - ((mkn^-k)^) 



(9) 



(10) 



in the logarithmic case. In the case of a 3D Coulomb gas, 
the interesting quantity will be e^/ LP' , which is indepen- 
dent of system size since (wk^-k) ^ L in that case. 



Standard Metropolis Monte Carlo simulations are car- 
ried out on the model Q at zero twist. An L x L square 
lattice with periodic boundary conditions is used and the 
system is kept electrically neutral at all times during the 
simulations. This is achieved by inserting dipoles with 
probability according to the Metropolis algorithm: An 
insertion of a negative or positive charge is attempted at 
random at a given lattice site, and an opposite charge 
is placed at one of the nearest neighbour sites to make 
the dipole. This is one move, accepted with probabil- 
ity exp{-AE/T) = exp[-(Hnow - Hoid)/^], and the 
sequence of trying this for all sites in the system once 
is defined as one sweep. If a charge is placed on top 
of an opposite one, the effect is to annihilate the exist- 
ing one. All simulations are performed going from high 
to low temperature and after simulating one system size 
L the sampled data are postprocessed using Ferrenberg- 
Swendsen reweighting techniques i2i 



A. 2D Coulomb gas 

We consider first the 2D Coulomb gas, which is known 
to suffer a metal-insulator transition via a Kosterlitz- 
Thouless phase transition. In this case, Monte Carlo data 
are obtained for L = 4 — 100 and for each L up to 200 
000 sweeps at each temperature is used. 

We start by taking the Hamiltonian 101 and comput- 
ing the mean square of the dipole moment, (s^), as a 
function of system size and temperature. A mean square 
dipole moment which is independent of system size in- 
dicates the existence of tightly bound dipoles and a di- 
electric or insulating phase. If the mean square dipole 
moment scales with system size, this demonstrates the 
existence of free unbound charges and hence a metallic 
phase. In other words, we expect in the low-temperature 
dielectric insulating phase no finite-size scaling of (s^), 
whereas we should expect (s^) oc L^C^) with a{T) < 2 at 
higher temperatures. Using an intuitive low density ar- 
gument, neglecting screening effectSj^ we can calculate 
the behaviour of (s^) to leading order in L, 



{Const. 
]^(T-Tkt)/T 



T<Tkt 

Tkt <T < 2Tkt 

2Tkt < T. 



ill) 



Hence, oi(T) is zero for low temperatures and a mono- 
tonically increasing function of temperature just above 
Tkt- Including screening effects in 2D shows that this 
conclusion still holds, however the temperature at which 
it occurs is determined by screening. 

Details of the simulations may be found in Ref. 
The result is shown in Fig. ^ where we have the mean 
square dipole moment for the 2D case both as a function 
of temperature for various system sizes, and as function 
of system size for various temperatures. From this we 



4 



may extract the scaling constant a{T) which is shown 
in the center panel of Fig. ^ A related method for us- 
ing dipole fluctuations to measure vortex-unbinding has 
recently been used in Refs. |2j 
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FIG. 1: The mean square dipole moment (s^) as a function 
of temperature (top panel), and system size (bottom panel) 
for the 2D Coulomb gas. The middle panel shows the scaling 
exponent a extracted from (s^) ~ L"'-"^^ 

Below a temperature T w 1.3, no scaling of (s^) is seen, 
consistent with a low-temperature dielectric phase. The 
temperature at which scaling stops is consistent with the 
known temperature at which the 2D Coulomb gas suffers 
a metal-insulator transition. 

Simulation results for the inverse dielectric constant 
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FIG. 2: Inverse dielectric constant taken at the smallest possi- 
ble wave vector in a finite system, k = (0, 2n/L), and plotted 
against temperature T for system sizes L = 10, 30, 50, 70 and 
100, for the 2D Coulomb gas. The decrease of towards 
zero becomes sharper with increasing L, consistent with the 
prediction of a discontinuous jump. Errorbars are given in the 
top and bottom curves, and omitted for clarity in the others. 



are shown for a selection of system sizes in Fig. |21 Since 
is expected to be discontinuous at in the limit 
k — ^ 0, we consider only the smallest possible wave vec- 
tor in each system, k = (0, 27r/L), and we see that the 
decrease of towards zero with increasing T indeed 
gets sharper as L grows. It is however difficult to decide 
from these plots alone whether or not the dielectric con- 
stant is discontinuous at Tc. The fourth order term in 
the expansion of the free energy, defined in eq. H1U|I . 
is therefore investigated in a corresponding manner and 
plotted in Fig. El 

We note that this quantity has a dip at a temperature 
which can be associated with the transition temperature. 
If this dip remains finite and negative as L approaches 
infinity, e^^ must exhibit a jump at Tc. The depth of the 
dip is shown in Fig. 2|for a variety of system sizes ranging 
from L = 4 to L = 100 and as a function oil/L. It clearly 
decreases with increasing L. However, from the positive 
curvature of the data in the log-log plot we may conclude 
that the depth remains nonzero when we extrapolate to 
l/L = 0, a conclusion reached by assuming power-law 
dependence of the depth on L. 

We can now subtract from the depth a constant cho- 
sen so as to linearize the curve in the log-log plot. This 
constant consequently corresponds to the depth when ex- 
trapolating the data to the thermodynamic limit l/L — 
0, and we find this to be 0.047 ± 0.005. 

By plotting the temperature at which the fourth or- 
der term has its minimum against l/L, we can follow a 
similar procedure as the above one. This is shown in 
Fig. |5l We linearize a log-log plot by subtracting a 
carefully chosen constant and end up with the number 
1.36 ± 0.04. This is nothing else than an estimate of the 
critical temperature of the 2D CG, and compares well to 
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FIG. 3: The coefficient 64 of the fourth order term of the 
expansion of the free energy, for the 2D Coulomb gas. The 
same systems are used in this plot as in Figure ^ and the 
depths decrease with increasing L. The important question is 
whether this dip vanishes at Tc or not. Errorbars are omitted 
but will be reintroduced in Figure |1| The oscillation at high 
T is due to noise from the reweighting. 
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FIG. 5: Temperature minimizing £4 as a function of inverse 
system size for the 2D CG. The values are plotted both on a 
linear scale and on a log-log scale (inset). This temperature 
reaches a nonzero value at L ^ 00 indicated by the positive 
curvature in the log-log plot. Extrapolation gives Tc = 1.36 ± 
0.04. 
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FIG. 4: Depth of the dip in the fourth order term shown in 
Figure|3for the 2D Coulomb gas. The data are obtained from 
simulations of system sizes ranging from L = 4 to L = 100 and 
plotted both on a linear scale (inset) and on a log-log scale. 
The positive curvature in the log-log plot clearly indicates a 
nonzero value of the depth when extrapolating to the limit 
L — > 00. 

earlier results.— The approach towards Tc is however a 
bit slow, making a precise determination of the critical 
temperature difRcult. This drawback was also noted by 
Minnhagen and Kim for the corresponding computations 
on the 2D XY model3 



B. 3D logarithmic system 

We may carry out the same type of analysis for the 
mean square dipole moment for a system of point charges 



interacting via a three-dimensional logarithmic bare pair 
potential (3D LG). For this system, much less is known. 
Such a system has recently been considered in the context 
of studying confinement-deconfinement phase transitions 
in the (2 + l)-dimcnsional abelian Higgs model. The 
results are shown in Fig. 

Qualitatively and quantitatively the results are the 
same in the 3D LG as for the 2D case. This strongly 
suggests that the 3D LG also has a low-temperature di- 
electric insulating phase separated by a phase transition 
from a high-temperature phase. In the low-temperature 
regime the charges of almost all dipoles are bound as 
tightly as possible, the separation of the charges corre- 
spond to the lattice constant. In the high-temperature 
regime the dipoles have started to separate, reflected by 
a scaling of (s^) ~ L"^"^^ with the system size. Since 
ct{T) = at low temperatures while a{T) 7^ in the high- 
temperature regime a non-analytic behaviour of a{T) is 
implied. This necessarily corresponds to a phase tran- 
sition in the vicinity of T « 0.3, a temperature which 
agrees well with Ref. ^3 where a critical value of Tc = 1 /3 
was obtained. 

Note that, although this simple type of analysis of the 
mean square dipole moment does not by itself suffice to 
determine the character of these phase transitions either 
in the case of 3D LG or 2D CG, it does suffice to shed 
light on the important issue of whether a low temperature 
insulating phase exists in the 3D LG as well. This is far 
from obvious, since the screening properties of a three- 
dimensional system of charges interacting logarithmically 
is quite different from that of a Coulomb system (in any 
dimension) It is therefore of considerable interest to 
repeat the analysis carried out for the 2D Coulomb gas to, 
if possible, determine the character of a metal-insulator 
transition in the 3D LG. 
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FIG. 6: Mean square dipole moment (s^) as a function of 
temperature (top panel) and system size (bottom panel) for 
the 3D system of point charges interacting with a logarithmic 
bare pair potential (3D LG). The middle panel shows the 
scaling exponent a extracted from (s^) ~ L"'^^ 




FIG. 7: Inverse dielectric constant taken at the smallest pos- 
sible wave vector in a finite system, k — (0,27r/L,0), and 
plotted against temperature T for system sizes L — 4, 10, 16, 
30, 40 and 56, for the 3D LG system. The decrease of 
towards zero becomes sharper with increasing L, consistent 
with the prediction of a discontinuous jump. However, the 
downward drift in the temperature at which the inverse di- 
electric constant starts decreasing rapidly is more pronounced 
than in the 2D CG case. Errorbars are given in the top and 
bottom curves, and omitted for clarity in the others. 
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T 

FIG. 8; The coefficient £4 of the fourth order term of the 
expansion of the free energy for the 3D LG model. The depths 
decrease with increasing L, and the important question is 
whether this dip vanishes at Tc or not. Errorbars are omitted 
but will be reintroduced in Figure |5| 



In Fig. dwe show the inverse dielectric constant for the 
3D LG as a function of temperature for various system 
sizes. It shows qualitatively the same behavior as for the 
2D CG in that the decrease of towards zero becomes 
sharper with increasing L. However, the downward drift 
in the temperature at which the inverse dielectric con- 
stant starts decreasing rapidly is more pronounced than 
in the 2D CG case. 

In Fig. |Slwe have plotted the fourth order coefficient 
against temperature for the 3D LG system, and the depth 



of the dip as a function of system size is shown in Fig. |51 
It would clearly have been desirable to be able to access 
larger system sizes than what we have been able to do in 
the 3D LG case, to bring out a potential positive curva- 
ture that was observed in the 2D CG case. From these 
results, it is unfortunately not possible to tell whether the 
depth of the dip remains finite and negative as L — > c» or 
if it vanishes. Hence, we are presently not able to firmly 
conclude that the inverse dielectric constant in the 3D 
LG experiences a discontinuity. 
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FIG. 9: Depth of the dip in the fourth order term shown in 
Figure |S| for the 3D LG. The data are obtained from simu- 
lations of system sizes ranging from L = 4 to L = 60 and 
plotted both on a linear scale (inset) and on a log-log scale. 
The lack of clear positive curvature in the log-log plot that 
was observed in 2D CG case makes the extrapolation to the 
limit L — > oo more difficult for the system sizes we have been 
able to access in 3D. 



The temperature locating the minimum in £4 as a func- 
tion of system size is shown in Fig. ^] for the 3D LG 
system. Extrapolation gives Tc = 0.30 ± 0.04. 
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FIG. 10: Temperature minimizing £4 as a function of inverse 
system size for the 3D LG system. The values are plotted 
both on a linear scale and on a log-log scale (inset). This 
temperature reaches a nonzero value at L ^ 00. Extrapola- 
tion gives Tc = 0.30 ± 0.04. 



exhibit quite different finite-size scaling of (s^) compared 
to the 2D CG caseiiiJ^iS^ The results are shown in Fig. 
im Note that the temperature dependence of the curves 
for all different system sizes are qualitatively different in 
the 3D CG compared to those in the 2D CG and the 
3D LG. This becomes particularly apparent upon con- 
sidering the i-dependence of (s^) for various tempera- 
tures, where the steepness of the curves increases with 
decreasing temperature, resulting in a scaling exponent 
a{T) (from (s^) ~ L"(^)) which decreases with increas- 
ing temperature. This is quite consistent with what is 
known for the 3D CG, namely that it exhibits a metal- 
lic state for all finite temperatures, equivalently it cor- 
responds to Polyakov's permanent confinement ^ii^ It is 
evident that the scaling results for (s^) for the 2D CG 
and the 3D LG are qualitatively and quantitatively the 
same, and that they are qualitatively different from those 
exhibited by the 3D CG. For low temperatures, (s^) seem 
to be increasing with temperature. This is only a vacuum 
effect, since vacuum configurations do not contribute to 
the measurement of (s^)i2i. This means that close to vac- 
uum, only configurations resulting from the insertion of 
one single dipole at the smallest possible distance will 
contribute. See also section ITVDI 

The inverse dielectric constant for the 3D CG is shown 
as a function of temperature in Figure 1121 with system 
sizes ranging up to L = 50. Here also, e^^ decreases 
from unity to zero, but the downward drift in the tem- 
perature at which e^^ deviates from unity seems to be 
even stronger than for the 3D LG model. Additionally, 
the decrease towards zero does not sharpen significantly 
with increasing L. 

We find a similar minimum in the fourth order term in 
the expansion of the free energy for the 3D CG, £4/^'^, 
shown in Fig. 1131 However, the dip vanishes as L ^ 00 
in the current model. This is clearly shown in Fig. 1141 in 
contrast to the Figs. 01 and El of the other two models. 

For completeness we have included in Fig. ^la plot of 
the temperature locating the minimum in £4 as a func- 
tion of system size also for the 3D CG. There is no phase 
transition to which this temperature is associated, and 
the stronger downward drift mentioned above is evident 
when contrasting this plot to Fig. EH of the 3D LG. 
The temperature is reduced by a factor 2 in the largest 
system considered in the 3D CG compared to the small- 
est whereas the variation is much smaller in the 3D LG. 
However, there is a weak curvature in the log-log ver- 
sion of Fig. 1151 Performing a similar extrapolation as we 
did for the other two models we end up with a "critical" 
temperature = 0.24 ± 0.04. 



C. 3D Coulomb gas 

In this subsection, we contrast the results of the 2D 
Coulomb gas and the 3D LG to those of the 3D Coulomb 
gas. The 3D CG is known to be in a metallic high- 
temperature phase for all finite temperatures and should 



D. Charge density 

Finally we present in Fig. ^1 the charge density for 
the three models considered. In all three cases the charge 
densities are independent of L and from these curves we 
can approximate the average separation Tmcan between 
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FIG. 11: Mean square dipole moment (s^) as a function of 
temperature (top panel) and system size (bottom panel) for 
the 3D Coulomb gas system of point charges interacting with 
a 1/r bare pair potential. The middle panel shows the scaling 
exponent a extracted from (s^) ~ L"'^'. 
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FIG. 12: Inverse dielectric constant taken at the smallest pos- 
sible wave vector in a finite system, k — (0,27r/L,0), and 
plotted against temperature T for system sizes L = 4, 8, 16, 
30 and 50, for the 3D CG system. The decrease of to- 
wards zero does not sharpen with increasing L, and there is 
a clear downward drift in the temperature at which devi- 
ates from unity. Errorbars are given in two of the curves, and 
omitted for clarity in the others. 
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T 

FIG. 13: The coefBcient £4 of the fourth order term of the ex- 
pansion of the free energy, for the 3D CG. The depths decrease 
with increasing L and seem to vanish as L ^ 00. Errorbars 
are shown for one of the systems for demonstration. 



the charges assuming uniform distribution, 



; Sum 



/V 



1/d 



(12) 



where d is the dimension. We consentrate on the (L- 
dependent) temperatures which minimize €4. In the two 
logarithmicaUy interacting models, rmean ranges from ^ 
4 for the smallest systems and up to ^ 8 for the largest. 
In the 3D Coulomb gas on the other hand, Tmcan remains 
close to L even for the largest system sizes meaning that 
the systems are close to their vacuum states at these 
temperatures. This strongly suggests that the features 



9 



o 

ft 

CO 

Q 



0.01 




0.001 



FIG. 14: Depth of the dip in the fourth order term shown 
in Figure 1131 for the 3D CG. The data are obtained from 
simulations of system sizes ranging from L = 4 to L = 50 and 
plotted both on a linear scale (inset) and on a log-log scale. 
It is clear that the dip vanishes in the thermodynamic limit. 
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system size for the 3D CG system. The values are plotted 
both on a linear scale and on a log- log scale (inset). 



we investigate are only extreme low-density effects in the 
3D CG model. Screening, which should take place at 
all temperatures in a system always being in a metallic 
state, is not possible in this limit. 

In the 2D CG and 3D LG models the situation is dif- 
ferent. The interesting temperature domains arc smaller 
and the charge densities are kept close to constant which 
in turn allows screening for the largest systems. 
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FIG. 16: Charge density Qsum/V^ plotted vs. temperature 
on log-log scales for the a) 2D CG, b) 3D LG and c) 3D CG 
models. The volume V corresponds to the total number of 
sites L"*. Note that Qsum/V^ is independent of system size L 
in all three cases. 



V. COMMENTS ON UNIVERSALITY 

In the 2D CG, the universal jump to zero of the inverse 
dielectric constant y^ is given bjii2iS& 



TT 



(13) 



Using the estimate for the critical temperature found in 
section IIV Al the value at Tc should, according to Eq. 
(IT2J), be = 0.86 ± 0.03. This is in agreement with 
Fig. 13 since it is in this region the curves seem to split. 

In Ref. [20I it was speculated that the finite negative 
value of the fourth order modulus (T4) « —0.130 could 



10 



be associated with a universal number. In the 2D CG, 
Vk ~ and ^ 2ti/L for large L, such that the 
modification A A/(A/27r) turns (jSJ into 



F{T) 



A2 



— 364. 



(14) 



This means that if e^^ corresponds to the helicity modu- 
lus (T), it is 3e4 that corresponds to the fourth order 
modulus (T4). It is interesting to notice that 864 = 
—0.141 ± 0.015 fits nicely with the value found in Ref. 
I20I speculated to be a universal number. One may there- 
fore speculate that the value of £4 at Tc is a universal 
number independent of Tc- Whether this is a sign of a 
true universality or a mere coincidence requires further 
investigation. 

One should also note that with this modification of A, 
the additional twist term in the XY-Hamiltonian Q be- 
comes •\/2Asin(27rj//L)/L. It seems natural to suggest 
that the net effect of a sine twist is given by its RMS- 
value, i.e. {l/L^l^s\T?{2'Ky/L)Ayf^ = l/y/2. This 
gives a net twist of A across the system, which is the 
same as in Ref. |20[ 

The universal jump of in the 3D LG is given by 
the flow equations derived in Ref. Q. In our units, this 
jump is predicted to be 



2 



(15) 



and by using the critical temperature found in section 
IIVBI this amounts to an in the interval (0.65,0.85). 
Since the different curves in Fig. Qdo not merge in the 
low-temperature regime, as they do in the 2D CG case, 
it is difficult to make a precise determination of the jump 
in the 3D LG based on these simulations. However, one 
can not rule out that the jump lies inside the interval 
mentioned. 



VI. 



CONCLUDING REMARKS 



In this paper we have considered various quantitities 
related to a possible phase transition in systems of point 
charges interacting with bare logarithmic pair potentials, 
in 2D and 3D. We have also carried out comparisons with 
the results obtained in the 3D Coulomb gas in some cases. 
The quantities we have focused on are the fluctuations 
of the dipole moment, (s'^), and the fourth order coeffi- 
cient of the free energy expanded in an appropriate twist. 
We have shown that the dipole moment fluctuations, as- 
sociated with the polarizability of the charge systems, 
has a scaling exponent a{T) defined by (s^) ~ L"(^) 
which is positive above some temperature and zero below 
this temperature for the 2D CG and the 3D LG cases, 
and is an increasing function of temperature. On the 
other hand, for the 3D CG case a{T) is finite positive for 
all temperatures we have considered, and is a decreasing 
function of temperature. This in itself strongly suggests 



that the 3D LG has statistical physics much more akin to 
the 2D CG than to the 3D CG. For the 2D CG we have 
demonstrated that the inverse dielectric constant expe- 
riences a discontinuous jump to zero at the phase tran- 
sition. This has been done by investigation of a series 
expansion of free energy using Monte Carlo simulations. 
The possibility of a universal value of the fourth order 
term proposed in Ref. [23 has also been commented on, 
and a possible agreement with this value has been ob- 
served. The method developed in this paper will apply 
to any gas of vortex loops or point charges with any in- 
teraction potential. We have applied it to the 3D LG. 
Although it would have been desirable to be able to ac- 
cess larger system sizes than what we have been able to 
in the present paper, the results we obtain for the 3D 
LG suggest that this model may also undergo a metal- 
insulator transition with a discontinuity in the inverse di- 
electric function at the critical point, in agreement with 
the renormalization group results of Ref. 
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APPENDIX A: DUALITY TRANSFORMATION 

The partition function for the XY-model with coupling 
constant J = 1 is 

^^n, /^e'3i:.™^(ve-2-T)^ (Ai) 

J 27r 

where the sum is over all links between lattice points, 
V6 = 9i — 6j and T(r) is the twist between the two 
lattice points sharing the link r. We will consider three 
spatial dimensions and comment on any differences in 
2D. Applying the Villain approximation, we get 

Z= f p0^e-^S^(^«-2.T-2.n)^^ (A2) 

n(r) is an integer- valued field taking care of the periodic- 
ity of the cosine. By a Hubbard-Stratonovich decoupling, 
one finds 

Z= f VeVv J2 e- (A3) 

The summation over n may now be evaluated using the 
Poisson summation formula. 



E 



(A4) 
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at each dual lattice point, yielding 



{1} 



2irilT-ilVe- 



2p' 



(A5) 



The field l(r) is integer- valued. Now, performing a par- 
tial summation on the second term in the exponent, 
the ^-integration may be carried out. This produces 
the constraint that 1 must be divergence-free, solved by 
the introduction of another integer-valued field such that 
1 = V X h. Note that h(r) is a scalar in 2D. The partition 
function is now 



{h} 



£,27ri(Vxh).T-5^(Vxh)2 



(A6) 



and we observe that h ^ h + Vi/) is a gauge transforma- 
tion. In two dimensions, the corresponding gauge trans- 
formation is h h + c, where c is a constant. Using 
Poisson's summation formula once more, we get 



{m} 



aEr2^i(Vxh)-T-^(Vxh)H27rih-m 



(A7) 



leaving h no longer integer- valued. The field m(r) is 
what corresponds to vortex excitations in the XY model. 
The gauge invariance of the theory produces the con- 
straint 't' ' ^) — for all configurations of m. 
Choosing for instance </) = V • m, it is clear that m must 
be divergence-free, i.e. the field lines are closed loops. In 
2D, the corresponding constraint is m = 0, indicat- 
ing an overall charge neutrality in the 2D Coulomb gas 
or zero total vorticity in the 2D XY model. 

By another partial summation, we are now left with 
a Maxwell term and a coupling term between the gauge 
field h and the current M(r) ee m V x T: 



{m} 



X,27rih-M-5L(Vxh) 



(A8) 



One may now perform a partial integration in the sec- 
ond term and use the gauge where V • h = 0, such that 
V X V X h = — V^h. Then, by going to Fourier space and 
completing squares, the h-integration becomes Gaussian. 
This leaves us with 



e— E 

{m} 



(A9) 



where V^e^"''" = e^"i''G'q and Zq is a constant. Defining 
the discrete Laplacian by 

AV(r) - E [/(^ + + /(r - e,) ~ 2/(r)] , (AlO) 



it is clear that Gq ~ ^2 ^rf — cos , denoting 

the number of space dimensions by d. Returning to real- 
space representation, we arrive at 



■E.,,., M(r,)V(|r.-r,|)M(r,) 



{m} 



the interaction being given by 



V{v) 



cosg^ 



(All) 



(A12) 



APPENDIX B: EXPANSION OF FREE ENERGY 



Consider the Hamiltonian 



= I'^rn^VijUij, 



(Bl) 



describing a 3D system of integer-valued currents m on 
a lattice interacting via the potential Vij = V{\ri — rj|). 
We impose periodic boundary conditions on the system. 
Perturbing the field m with a transversal twist turns ljBl|) 
into 



H= i^(m + V X T)^ViJ{m + \/ x T)j 



(B2) 



We let the linear system size be L and define the discrete 
Fourier transform by 



/q = ^/(r)e'^ 



(B3) 



where r = [ux, ny,nz) and nt ~ 0, i — 1. The inverse 
transform is 



(B4) 



where q = ^{k^^ky, k^) and k, = -L/2 + 1, L/2. N 
is the number of lattice sites. Let us also define by 
^i^giiq r ^ giiq rg^^^ where is a lattice derivative. 
In Fourier representation the Hamiltonian becomes 



12 



2iV 



(B5) 



For later use, we calculate the derivative of H, which is 



BM 1 
qi 



(B6) 



We also note that 



e'"'"£'^'"'Q-q,Q^,Kjx'5q,+q.,0 



(B7) 



is independent of m and that all higher order derivatives are zero. 

The free energy is given hj F = —T In Z, where the partition function is 



{m} 



summing over all possible configurations of m. By Taylor expansion of the free energy in the twist, we get 



(B8) 



F{T)-F{Q) = Y,Y. 



OF 



5T«(ri) 



T=0 



d^F 



a, (3 ri ,r2 



5T«(ri)aT/3(r2) 



r"(ri)r'3(r2) + ... (B9) 



T=0 



Note that F(T = 0) refers to the free energy of the unperturbed system described by Hq. By writing each term in 
the series in Fourier representation, one finds the equivalent expansion in Fourier components of the twist, i.e. 



F(T)-F(0) = ^^|^ ^^q.+EE 



qi T=0 



a,l3 qiq2 '^-'qi^-'q2 



qi q2 



+ . 



T=0 



(BIO) 



The first derivative becomes 



dF j_ ^ dl£_ _H/T ^ / dH 



(Bll) 



Proceeding, we find 



d'^F 



1 / dH \ dH 



1 / dH dH 

OT-dTi ~ T \dfl/ \a7f / " OT-dTg ' T \dT- dT^ 
for the second derivative and 

d^F 



1 



dH \ dH 



dH 



dT^J \dTi / \dT, 
dH \ I dH dH 



dH 



/ dH dH dH 

^[dTidrgdrl 

dH \ / dH dH 



dT^JTl 
' dH dH 



dTi \dT-^dTl 



dTl/ XdT^^dTi 



(B12) 



(B13) 
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for the third. We have exploited the fact that third derivatives of H vanishes. The fourth derivative is found to be 



dT-^dTgdTldTl 




dH \ 


/ dH 


dH 


dTl) 


\dK 


dT^ 


dH \ 


1 dH 


dH 


dTlJ 






dH \ 


1 dH 


dH 




\9TS. 


dTi 



Remembering that 



/ dH dH 
\dfldfl 



dH 



dH dH dH ' 

djf^dJZdfl, 

dH dH dH ' 
dT^^OTgdTl^ 

dH dH \/ dH dH ^ 
dfldfZ/\dTtdfl^ 

dH dH dH dH 
or- dTl dTl dTl 



(B14) 



dT- 



T=0 



(B15) 



it is straightforward to write the derivatives at zero twist as m-correlators. However, in many cases these expressions 
may be simplified further. If the sum over all possible configurations {m} is symmetric around zero, one finds that 
all odd-order correlators are zero, resulting in 



dF 



dT^ 



d^F 



T=0 



dT-^dTi,dTl 



0. 



(B16) 



T=0 



Furthermore, since Vij = V{\ri 
subject to relations like 



i.e. we have a translationally invariant system, the even-order correlators are 



and 



Thus, we find 



+q2+q3+q4,o- 



d^F 



e''""£'^^^Q^q,Q^q.V'q,<5q,+q„0 



&T^,dT^. 

for the second derivative and 

d-^F 



T=0 



N 



("l-qi"^-q2)) 



NT 



(B17) 
(B18) 

(B19) 



dT»dTidTldTl 



^^aa^.pP^.r^^X^SQa_^^ qP_^^ Qr_^^ QV_^^ y^^ y^^ y^^ y 



"•it 



T=0 



N^T^ 



+ ( 
+ ( 



u A 

"^-qi"i-q4 



l) '5qi+q2,0 ^q3+q4,0 



) ^qi+q3,0 0q2+q4,0 
) '^qi+q4,0 '5q2+q3,0 
) ^qi+q2+q3+q4,o} 



(B20) 
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for the fourth. These expressions may also be apphed to a gas of point charges in 2D or 3D, that is when m is a scalar 
field. One way to do this is by replacing V x T in (|ij2l) with its z-component e^'^'^'A'^T'^, with the consequence that 
the greek letter summations may be taken over x and y only. For the second derivative, this results in 



qi '^qi"qi+q2,0 



T=o 



N 



NT 



(B21) 



where we have applied Vq — V-q. We recognize the 
paranthesis as the Fourier transform of the inverse dielec- 
tric response function e~^(qi) in the low density limit. 
Note that the factor 

^..a^.,/,Q.^^QP^ = g-^g-q^ (^i ^1^) (B22) 

is a projection operator times Qq^^Q-q^ ^ qfx + Qiyi re- 
flecting the transversality of the twist. 

To arrive at eq. we chose the twist (Ujl and com- 
puted the sums appearing in the expansion l|B10|l for 
both the second and fourth order term. The sum over di- 
rection is trivial, since our twist points in the x-direction. 
The sum over momenta is also managable, since has 
nonzero values only for q = (0,±27r/L). This sum gives 
two contributions in the second order term, due to the 
restriction (Jqj+qj^o- The same argument results in four 
contributions for the three terms in (|B20I) being a prod- 
uct of two second order correlators. The term containing 
a fourth order correlator will give six contributions due 



to the restriction (Jq^+q^+qj+q^^o- 

APPENDIX C: HIGHER ORDER TERMS 

Using the method described in this paper involves ex- 
trapolation to i ^ cxD and deciding whether or not the 
fourth order term in the expansion (|B10|) goes to zero or 
to a finite nonzero value. This procedure could in some 
cases be difficult. However, if the fourth order term had 
turned out to be zero in the thermodynamic limit, it 
would not necessarily mean that the second order term, 
the inverse dielectric response function, would have to go 
continuously to zero. In fact, if one were able to prove 
that the fourth order term is negative or zero, one could 
go on to investigate the sixth order term instead. If it 
then turned out that the value of the sixth order term 
was hard to establish, one could in principle repeat the 
procedure and go to higher order terms. We therefore in- 
clude the sixth derivative here. To simplify calculations, 
we work with a twist in the x-direction only. 
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Note that we are allowed to permute the momenta qi, qe, since these are summed over in the free energy expansion. 
Assuming vanishing odd-order correlators and imposing that m is a scalar field gives 



ri'Tx fj'Tx rl^x rl'Tx fl'Tx Pi'Tx 

qi q2 q3 q4 q5 q6 



X <; 12 (m_q^m„q2 



NT 



-qs '"'-q4/ 



'^qi+q2,0 '5q3+q4,0 (^qs+qg.O 
(TO-q^m-q^m-qgrn-q^m-q^m-qJ (5qi+q2+q3+q4+q5+qg,0 (C2) 

-q4'™-q5''^-q6) '^q3+q4+q5+q6,0 
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